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FORTRAN PROGRAM FOR QUASI -THREE -DIMENSIONAL CALCULATION 


OF SURFACE VELOCITIES AND CHOKING FLOW 
FOR TURBOMACHINE BLADE ROWS 
by Theodore Katsanis 
Lewis Research Center 

SUMMARY 

A quasi -three -dimensional compressible flow analysis is reported in NASA TM X- 
1394 for axial flow turbines. This analysis has been generalized herein to allow for 
mixed or radial flow and for nonuniform inlet conditions. The velocity gradient method 
is used, with one velocity gradient equation used for radial equilibrium and another ve- 
locity gradient equation for blade -to -blade variation. In addition to a description of the 
quasi -three -dimensional flow analysis, a FORTRAN IV computer program is included 
for the solution of these equations. 


INTRODUCTION 

The aerodynamic design of turbomachinery blades requires the determination of the 
blade surface velocity distribution. Also, it is often necessary to know the choking mass 
flow. In many blade designs there are significant velocity gradients both from blade to 
blade and from hub to tip; this condition necessitates consideration of three-dimensional 
effects. 

One of the useful techniques for calculating surface velocities where three- 
dimensional effects are of importance is the velocity gradient (stream filament) method. 
The general velocity gradient equation (ref. 5) determines the velocity variation in any 
direction. In particular, the velocity gradient equation can be reduced to special cases 
to determine both the blade to blade and hub to tip variation in velocity. A combination 
of the velocity variation in two directions with a specified mass flow will determine the 
velocities at a passage cross section. This method works well in a well-guided passage. 
Reference 1 applies this method to axial flow turbines. The computer program (CTTD) 



in reference 1 has been used at Lewis Research Center for 17 years and has long been 
available to industry for their own use in turbine design. 

A new program, CHANEL, has been written. CHANEL is more general than CTTD 
and can obtain quasi -three -dimensional solutions in any well-guided channel. Some of 
the things that can be handled by the CHANEL program that could not be handled before 
are nonuniform inlet temperature, pressure, and prewhirl, and nonaxial flow where 
meridional flow angle, meridional streamline curvature, and radius can vary as desired 
from the hub to tip. Also, output has been clarified, and additional output information is 
given. 

This report gives a description of the analysis procedure and of the use of the 
CHANEL program. The CHANEL program listing is also given. 
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SYMBOLS 


coefficient, eqs. (Bl) and (B2) 
coefficient, eqs. (A3) and (A4) 
coefficient, eqs. (Al) and (A2) 
coefficient, eqs. (Bl) and (B2) 
coefficient, eqs. (A3) and (A4) 
coefficient, eqs. (Al) and (A2) 
coefficient, eqs. (Bl) and (B2) 
coefficient, eqs. (Al) and (A2) 

specific heat at constant pressure, j/(kg-K) (ft-lbf/(slug-°R)) 

enthalpy, j/kg ((ft-lbf)/slug) 

meridional streamline distance, meters (ft) 

distance from hub along normal to meridional streamline, meters (ft) 


distance from suction surface along normal to blade -to -blade streamline, 
meters (ft) 

pressure, kg/meter^ (lbf/ft^) 
distance along an arbitrary curve 
gas constant, j/(kg-K) (ft-lbf/(slug-°R) 
radius, meters (ft) 
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Subscripts: 


radius of curvature of meridional streamline, meter (ft) 

radius of curvature of blade -to -blade streamline, meter (ft) 

radius at intersection of normal to meridional streamlines with hub (tip), 
meters (ft) 

temperature, K (°R) 

tangential component of absolute velocity, meters/sec (ft/sec) 

velocity relative to blade, meters/sec (ft/sec) 

critical velocity relative to blade, meters/sec (ft/sec) 

meridional component of velocity relative to blade, meters/sec (ft/sec) 

velocity relative to blade at midchannel between suction and pressure sur- 
faces and at the mean section between hub and tip, meters/sec (ft/sec) 

tangential component of velocity relative to blade, meters/sec (ft/sec) 

mass flow through the channel, kg/sec (slugs/sec) 

axial coordinate, meter (ft) 

angle between meridional streamline and axis of rotation, rad, see fig. 3 

angle between relative velocity vector and meridional plane, rad, see 
fig. 3 

specific heat ratio 

relative angular coordinate, rad, see fig. 3 

n o 

prerotation (rVg)^, meter /sec (ft /sec) 

3 3 

density, kg/meter (slugs/ft ) 

rotational speed, rad/sec 


i inlet or upstream 

isen isentropic 

loss difference between isentropic and actual 

Superscripts: 

f absolute stagnation condition 

,? relative stagnation condition 
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METHOD AND ASSUMPTIONS 


The objective of the analysis method is to calculate the quasi -three -dimensional ve- 
locity distribution satisfying continuity at a given channel orthogonal surface (see fig. 1). 
The mass flow may be specified for the calculation, or the maximum (choking) mass 
flow for that channel orthogonal surface may be calculated. The velocity variation on the 
orthogonal surface is given by differential equations for the rate of change of velocity. 
These equations are called velocity gradient equations. One velocity gradient equation 
gives the blade -to -blade variation, and the other gives the hub-to-tip variation. The ve- 
locity gradient equations are solved simultaneously with the condition of either a spec- 
ified or choking weight flow. This determines the blade surface velocities on a particular 
channel orthogonal surface. When this is done for several orthogonal surfaces, a ve- 
locity distribution over the blade surface is determined. 

The velocity distribution can be obtained in this way for the guided channel formed 
by the portion of the passage where the blade -to -blade orthogonals extend from suction 
to pressure surface. The guided channel will not cover the entire suction surface. To 
obtain the velocity distribution on the uncovered portion of the blade, the location of the 
stagnation streamline would have to be known. Therefore, CHANEL cannot be used to 
obtain velocities on the uncovered portion of the blade. (However, other methods are 
available for this problem, e. g. , refs. 3 and 4.) 

The basic simplifying assumptions used in deriving the equations are the following: 

(1) The flow relative to the blade is steady. 

(2) The fluid is a perfect gas with constant c . 

(3) The fluid is a nonviscous gas. 

(4) The midchannel line is a streamline, hereinafter referred to as the midchannel 
streamline. 

(5) There is a linear variation of streamline curvature along an orthogonal or there 
is a linear variation of radius of curvat ure along an orthogonal. An option is provided 
for this assumption in the program. 

(6) There is no change in radius along a blade -to -blade orthogonal. (The change in 
radius usually has a negligible effect on the solution.) 

(7) There is no change in flow angle along a blade -to -blade orthogonal. (The change 
in flow angle usually has a negligible effect on the solution. ) 

The method used here is very similar to that in reference 1. The main difference is 
that the assumptions of axial flow and uniform inlet conditions have been dropped. The 
analytical equations and the details of the solution procedure are given in appendix A. 
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DESCRIPTION OF INPUT AND OUTPUT 


The computer program requires, as input, operating conditions, gas constants, and 
flow passage geometry. The usual procedure for obtaining this input includes the follow- 
ing steps: 

(1) Define inlet absolute stagnation temperature and density, prerotation from hub 
to tip, mass flow per passage, rotational speed, and gas constants. 

(2) Lay out graphically (or determine analytically) the meridional profile of the pas- 
sage and determine hub to tip meridional normals as desired. Each meridional normal 
is a separate case requiring a separate set of input. Any number of cases may be in- 
cluded in a single computer run. From the meridional profile, the distance along the 
meridional normal from the hub, the radius from axis of rotation, the meridional flow 
angle, and the meridional streamline curvature must be obtained several points from hub 
to tip. 

(3) Lay out the blade -to -blade passage at hub, mean, and tip to determine blade -to- 
blade flow angle and orthogonal distance, and the suction and pressure surface curva- 
tures. 

(4) Prepare input sheets, one input sheet for each case. 


Input 

The input data sheet is shown in figure 2. The quantities shown are for a sample 
problem. The output for this sample problem is presented later. Certain items (marked 
*) are discussed in greater detail in the section Instructions for Preparing Input. 

The input variables are as follows: 

JX* controls blade -to -blade curvature variation assumption 

JX = 1, linear variation of curvature 
JX = 2, linear variation of radius of curvature 
JZ* controls type of solution obtained 

JZ = 1, subsonic solution 
JZ = 2, supersonic solution 
JZ = 3, choked flow solution 
KR1* controls use of card starting with GAM 
KR1 = 0, omit card with GAM 
KR1 = 1, supply card with GAM 
KR2* controls use of cards starting with RHOIP 
KR2 = 0, omit cards with RHOIP 
KR2 = 1 , supply one card with RHOIP 
KR2 = 3, supply three cards with RHOIP 
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NSP* 


number of meridional input points between hub and tip for the NAMELIST 
arrays (see ref. 8) (minimum of 2, maximum of 25) 

GAM specific heat ratio, y 

AR gas constant, j/kg-K (ft-lb/slug-°R) 

OMEGA rotational speed, tu, rad/sec (Note that w is negative if rotation is in the 
opposite direction of that shown in fig. 3. ) 

WTFL mass flow per blade, kg/sec (slug/sec) (not required for choked solution) 

3 3 

RHOIP inlet stagnation density, p!, kg/meter (slug/ft ) 

PLOSS* fractional loss of relative stagnation pressure (1 - (p"/p. M )) 

IS Gil 

NBB* length of streamline normal from blade to blade, meter (ft) (fig. 6) 

CS* blade surface curvature at intersection of orthogonal with suction surface, 

1/meter (l/ft) (fig. 6) (Sign is negative if surface is concave downward.) 
CPR* blade surface curvature at intersection of orthogonal with pressure surface, 
1/meter (l/ft) (fig. 6) (Sign is negative if surface is concave downward.) 

The remaining input is given in NAMELIST format as shown on the last two lines of 
figure 2. These variables are all arrays, with NSP items in each array. NSP must be 
at least 2 but not over 25. All the NAMELIST arrays are initialized to zero by the pro- 
gram so that zero input items may be omitted. Any input given in NAMELIST format 
remains unchanged from the previous case unless given again as input. 

NMERID* array of distances from hub along normal to the meridional streamlines, 

meter (ft) (These values must be in increasing order. The first value must 
be zero, and the last value must be the total distance along the normal from 
hub to tip. A meridional streamline is the projection on a plane through the 
axis of rotation of a streamline midway between the two blade surfaces.) 
(See fig. 5. ) 

TEMPIP* array of inlet stagnation temperatures T! for streamline corresponding to 
the NMERID array, K (°R) 

LAMBDA* array of values of prerotation A for streamlines corresponding to the 

2 2 

NMERID values, meter /sec (ft /sec) 

RADIUS* array of radii r corresponding to the NMERID array, meter (ft) (fig. 5) 
CURV* array of meridional streamline curvatures l/r c corresponding to the 
NMERID array, l/meter (l/ft) (fig. 5) 

ALPHA* array of meridional streamline angle a corresponding to the NMERID 
array, deg (fig. 5) 

BETA* array of relative flow angle /3 measured from the meridional plane, positive 
in direction of rotation, midchannel only, and corresponding to the NMERID 
array, deg (fig. 6) 
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Instructions for Preparing Input (see fig. 2) 


Units of measurement. - The International System of Units (ref. 2) is used through- 
out this report. However, the program does not use any constants which depend on the 
system of units being used. Therefore, any consistent set of units may be used in pre- 
paring input for the program. In particular, U.S. Customary Units may be used, and 
these are specifically noted in the list of input variables. 

Title cards . - The first card is a general title card for the entire run. Several 
cases may be submitted on a single run. And each of these cases may be labeled by 
means of the label card. All 80 columns may be used. 

Option code card. - The next card provides for several options in running the pro- 
gram. The first option (JX) is to choose an assumption of how the blade -to -blade stream- 
line curvature varies. The usual assumption is linear variation of curvature (JX = 1). 
This assumption must be used if the blade curvature is zero on either blade or if the 
curvature changes sign between the blades, since the curvature must be zero (r c = °o) at 
some point. On the other hand, there are cases when the assumption of linear variation 
of curvature is not satisfactory. For example, if the curvature is very large on one 
surface, and very small on the other, it may be much more reasonable to assume linear 
variation of radius of curvature (JX = 2). Ordinarily, there is not a large change in 
curvature, and in this case either assumption will give good results. 

The next option (JZ) specifies the type of solution desired. If there is a solution for 
a given mass flow through the channel, there are usually two distinct solutions. These 
two solutions have two different velocity levels. The solution at the lower level has 
mostly subsonic velocities, although some velocities may be supersonic. On the other 
hand, the solution at the higher level has mostly supersonic velocities although there may 
be some subsonic velocities. The solution with lower velocities will be referred to as 
the subsonic solution and the greater solution as the supersonic solution. The option 
JZ = 1 (=2) will result in obtaining the subsonic (supersonic) solution. If the maximum 
possible or choking mass flow solution is desired, use JZ = 3. The program will com- 
pute the choking mass flow and the corresponding velocities. If the mass flow specified 
as input is larger than the choking mass flow, the result will be the same as when JZ = 3. 

For the first case of a given run KR1 must be 1. Then each item on the card start- 
ing with GAM must be specified. For successive cases, if nothing on the GAM card 
changes from the previous case the card is omitted and KR1 is 0. If it is desired to 
change some items, then let KR1 = 1 and completely fill in the GAM card. 

The next option (KR2) specifies whether information is given on the RHOIP cards. 

On the first case of a given run, KR2 must be either 1 or 3 corresponding to 1 or 3 
RHOIP cards. For successive cases, if nothing on the RHOIP cards changes from the 
previous case, the cards are omitted and KR2 is 0. Use one RHOIP card if all quantities 
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are constant from hub to tip. When 3 RHOIP cards are needed, the quantities on these 
cards are given at hub, mean, and tip, in that order, where the mean is taken to be mid- 
way along the meridional normal between hub and tip. 

The next quantity (NSP) specifies the number of meridional input points for the 
NAMELIST arrays. A great deal of flexibility is provided here. NSP must always be at 
least 2, since information must be provided at hub and tip. Any additional points, up to 
a total of 25, are entirely at the discretion of the user, both as to the number and the lo- 
cation. 

Loss assu mption. - A loss correction may be made by specifying the fraction, 
PLOSS, of the relative stagnation pressure which has been lost. This is equivalent in a 
simplified calculation like this to specifying the fraction of the passage which is occupied 
by the displacement thickness of the boundary layer. PLOSS is zero if there is no loss 
correction. 

Meridional p rofile layout. - The meridional profile should include hub, shroud, and 
mean line, as well as blade leading and trailing edge, as shown in figure 5. Any number 
of meridional normals can then be drawn. The meridional normals should be normal to 
the hub, mean, and tip lines. If a meridional solution has been previously obtained, the 
meridional normal should be normal to the meridional streamlines. Each orthogonal will 
be a separate case using one input sheet. The quantities which are obtained from the 
meridional layout are NMERID, RADIUS, CURV, and ALPHA. Normally these four 
quantities are obtained at hub, mean, and tip. However, if a meridional .solution has 
been previously obtained more values can be given. 

NMERID is the distance along the meridional normal, so the first value must be 
zero and the last value must be the total length of the meridional normal. Any desired 
intermediate points may be added, up to a maximum total of 25. The total number of 
points is NSP. After the desired points have been chosen, all the lengths are measured 
from the hub. Then the corresponding radii are measured for the RADIUS array. The 
angles from horizontal, in degrees, go in the ALPHA array. If the hub and tip are 
straight, the wall curvature (CURV) is zero, so that no entry need be made for CURV. 

The example input in figure 2 is for straight axial flow so that both CURV and ALPHA 
are zero, and hence omitted from the input. If the hub or tip are not straight, the wall 
curvature should be measured at the ends of the meridional normal. It is suggested that 
the meridional streamline curvature be assumed to vary linearly between hub and tip, 
unless a meridional streamline solution is available. The curvature is considered to be 
positive if the streamline is concave upwards. 

After the NMERID distances have been determined, the inlet conditions of inlet stag- 
nation temperature (TEMPIP) and prerotation (LAMBDA) can be specified at corres- 
ponding points. TEMPIP and LAMBDA must each have NSP values given, although these 
values may be constant. 
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Blade -to -bla de channel layouts . - The blade -to-blade channel layouts are usually the 
most difficult step in obtaining the input data. If the hub and tip are straight (in the me- 
ridional profile layout), the blade -to -blade layouts can be laid out as a flat surface. If 
either the hub or tip is not straight, then the corresponding blade -to -blade surface is not 
developable and cannot be laid out as a flat surface with correct angles and linear meas- 
urements. There are several ways that a nondevelopable blade -to -blade surface of revo- 
lution can be laid out. However, these layouts are distorted, so that corrections must 
be made to determine orthogonal directions; also, distance measurements from the lay- 
out may require a graphical integration procedure. With the proper corrections, the de- 
sired blade -to -blade normals can be determined, and the lengths calculated with sufficient 
accuracy. Since this portion of the graphical procedure can be done in any of several 
ways, no further explanation is given here for determining blade -to -blade normals when 
hub or tip are not straight. 

Blade -to -blade layouts at hub, mean, and tip should be drawn as shown in figure 6. 
The layouts in figure 6 correspond to the blade shown in figure 5. The intersection point 
of the meridional normal with the mean line on the meridional layout (marked (T) in fig. 5) 
should be marked on the midchannel line on the mean blade -to -blade layout (marked (T) 
in fig. 6). This determines the reference 9 coordinate for the meridional normal. This 
reference 6 coordinate can be marked on the hub and tip layouts as shown in figure 6. 
Now the intersection point of the meridional normal with the tip on the meridional layout 
is marked on the tip blade -to -blade layout (marked (T) in figs. 5 and 6). The correspond- 
ing hub point is marked (J3)in figures 5 and 6. Now the true blade -to-blade normals are 
drawn to be normal to suction and pressure surfaces and to pass through the points (T) , 
(IT) , and (T) on the hub, mean, and tip layouts as shown in figure 6. Then the blade-to- 
blade normal length (NBB) is measured or calculated from each layout. Also, the blade 
suction surface curvature (CS) and blade pressure surface curvature (CPR) are meas- 
ured at the ends of the blade -to -blade normal. Blade surface curvatures are important, 
but they are difficult to measure accurately; therefore, extra care should be taken to 
measure blade curvatures accurately. Finally, the angle (BETA) of the mid channel line 
from axial is measured (with correction if necessary) at the hub, mean, and tip. The 
angle BETA is positive if the tangential component of velocity is in the direction of ro- 
tation. Since BETA is given as NAMELIST input, values are usually given at hub, mean, 
and tip, but more values must be given if NSP is greater than 3. 


Output 


An example of the output obtained is given in figure 7. 
the output given on the input form shown in figure 2. U.S. 


This output corresponds to 
Customary Units of measure - 
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ment are used. Each case will normally be printed on one page. The first output is a 
listing of all the input, pretty much in the same format as on the input sheet (fig. 2). All 
input items are given, including those which are repeated from the previous case. 

After the input listing, the velocities across the passage from blade to blade and from 
hub to tip are printed. Then the critical velocity ratio W/W cf across the passage is 
printed. If the passage is choked, the choking mass flow is also printed. 

Finally, information on the meridional streamline spacing is printed. Information is 
given at 10 equally spaced points from hub to tip. The information is given in three ways. 
First is the integrated mass flow fraction between the hub and the given point. The next 
column gives the ratio of the actual streamline spacing to uniform spacing. Finally, the 
streamline spacing (or normal stream sheet thickness) for 1 percent of the mass flow is 
given. This last item is useful if it is desired to obtain a blade -to -blade solution using 
the TSONIC or TURBLE program of reference 3 or 4. 


Error Messages 

A number of error messages have been incorporated into the program. These error 
messages are listed here. Where necessary, suggestions for finding and correcting the 
error are given. 

NSP MUST BE GREATER THAN OR EQUAL TO 2. 

NSP is the number of points in the NAMELIST arrays. Since the first point must be 
at the hub and the last point at the tip, NSP must be at least 2. 

NSP CANNOT BE GREATER THAN 25. 

JX MUST BE EITHER 1 OR 2. 

WHEN JX = 2, CS AND SPR MUST HAVE THE SAME SIGN. 

When JX = 2 it is assumed that the radius of curvature of blade -to -blade streamlines 
varies linearly. This is not possible if CS and CPR are of opposite sign, or if either CS 
or CPR is zero. 

A SOLUTION CANNOT BE OBTAINED AT THIS STATION. 

This message is printed when the solution procedure breaks down and no solution can 
be found. The input as listed should be examined carefully for errors. 

SPLINT USED FOR EXTRAPOLATION. EXTRAPOLATED VALUE = X.XXX 

This message is followed by a printout of the SPLINT input and output arguments. 
SPLINT is used for interpolation of the input values in the NAMELIST arrays. If 
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SPLING is used for extrapolation it must be due to an input error, most likely in the 
NMERBD array. 


CHANEL COMPUTER PROGRAM 

The program consists of the main program and the three subroutines SPLINT, PABC, 

and CONTIN. Subroutine SPLINT is used for interpolation by a spline fit curve, and is 

described more fully in reference 5 and 6. Subroutine PABC calculates the coefficients 

o 

A, B, and C of the parabola y=Ax +Bx + C passing through three given points. Subroutine 
PABC is called by both the main program and subroutine CONTIN. Subroutine CONTIN 
finds the value of W m ^ mean which will give either the desired value of mass flow w 
or the choking mass flow, depending on the input value of JZ. 


Main Program 


COMMON SRW 

DIMENSION NMERID(25) , T EM P I P ( 25 ) , L AMBDA ( 2 5 > , RAD I US ( 25 ) , CURV ( 2 5 ) , 

1 ALPHAI25) ,BETA(25 ) ,AAA( 25) 

DI MENS I ON NMER(9),TIP(9),LBDA(9),RAD(9),CRV(9),ALPH(9),BTA(9), 

1 AA<9),BB(9),CC(9) , W ( 9 ) , WB TB ( 9 ) 

DIMENSION RH0IPI3) , P LOSS ( 3 ) , NBB ( 3 ) ,CS(3),CPR(3) , SUM ( 3 ) , CPTI P ( 3 ) , 
I WCR(3),X(3),Y(3),BBB( 3) 

DIMENSION WWCR(3,3),WFIN(3,3) 

DIMENSION AABI9.3) 

REAL NMER I D, LAMBDA, NMER, LBDA, NBB 
INTFGER BTB , HMT , SRW 

NAMELIST /N AM 1/ NMER I D, TEMPI P, LAMBDA, RADIUS, CURV, ALPHA, BET A 
DO 5 1=1,200 
5 NME R I D ( I ) =0. 

WRirE(6,1000) 

READ (5,1010) 

WR ITE(6«1010) 

10 READ (5,1010) 

WRITE( 6, 1010 ) 

READ (5,1020) JX , J Z , K R 1 , KR2 , NSP 
WRITE (6,1100) 

WRITE(6,1020) JX»JZ,KR1,KR2,NSP 
IF(NSP.GT.l) GO TO 20 
WR I TE ( 6 , 1 500 ) 

STOP 

20 I F ( NSP . LE . 25 ) GO TO 30 
WRITF( 6, 1510 ) 

STOP 

30 WRITE (6,1110) 

IF(KRl.EQ.l) READ (5,1030) GAM , AR , OMEGA , WTFL 
WRITE* 6, 1040 ) GAM, AR, OMEGA, WTFL 
WRITE(6,1120) 

I F ( KR2 • GT .0 ) READ (5,1030) RHO I P ( 1 ) , PLOSS (1) , NBB (1) , CS ( 1 ) , CPR ( 1 ) 
I F ( KR2 - 1 ) 60,40,35 
35 CONTINUE 

READ (5,1030) RHO I P ( 2 ) , P LOSS ( 2 ) , NBB ( 2 ) , C S ( 2 ) , CPR ( 2 ) 
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o o o o n n 


REAU (5,1030) RHO I P ( 3 ) , P LOS S ( 3 ) , NB8 ( 3 ) , CS ( 3 ) ,CPR(3 ) 

GO TO 60 

40 DO 50 HMT = 2 * 3 

RHOIP(HMT) = RHO I P ( HMT- 1 ) 

PLOSS(HMT) = PLOSS(HMT-l) 

NBB(HMT) = NBB(HMT~1 ) 

CS(HMT) = CS(HMT-l) 

50 CPR (HMT ) = CPR(HMT-1 ) 

60 WRITE (6, 1130) ( RHO IP ( HMT ) , PL OSS ( HMT ) , NBB ( HM T ) , C S ( HMT ) , C PR ( HMT > t 
1 HMT = 1 » 3 ) 

READ ( 5 f N AMI ) 

70 WRITE ( 6 * 1140) 

WR I TE ( 6 t 1 040 ) (NMER ID (I), 1=1, NSP) 

WRITE! 6, 1 150) 

WR I TE ( 6 , 1040) ( TEMPI P( I ) ,1=1 , NSP) 

WR I T E ( 6,1160) 

WR I TE ( 6,1040) ( L AMBDA ( I ) ,1=1, NSP) 

WRI TE I 6 , 1 1 70 ) 

WRI Tfc( 6, 1040) (RADIUS ( I ) ,1 = 1, NSP ) 

WRITE( 6*1180) 

WRI TE ( 6, 1040 ) ( CUR V (1), 1=1, NSP) 

WR I TE ( 6 , 1 1 90 ) 

WR I TE ( 6 , 1040 ) (ALPHA (I), 1 = 1, NSP) 

WRI T E ( 6,1200) 

WR I TE ( 6 , 1040 ) (BETA (I), 1=1, NSP) 

IF (JX.EQ.l) GO TO 78 
IF (JX.EQ.2) GO TU 72 
WRITE (6,1000) 

GO TO 10 
72 DO 74 1=1,3 

74 IF (CS ( I )*CPR( I ) • LE. 0 • ) GO TO 76 
GO TO 78 

76 WRITE (6,1530) 

WRITE (6,1000) 

GO TO 10 

ALL INPUT HAS BEEN READ 

CALCULATE CONSTANTS AND INTERPOLATE MERIDIONAL VARIABLES AT 
EIGHT EQUAL INTERVALS FROM HUB TO TIP 

78 CP = AR/ ( GAM-1 * ) *G AM 
EXPON = 1. /(GAM-1.) 

TGROG = 2. *GAM*AR/(GAM+1 . ) 

DO 80 1=1,9 

80 NMER(I) = FLOAT ( I-l)*NMERID(N$P)/8. 

DELTA = N MER I D ( NSP ) / 8 . 

NINE = 9 

CALL SPLINT ( NMERID, TEMPI P, NSP, NMER, NINE, TI P,AAA) 

CALL SPLINT ( NMER ID, LAMBDA, NSP, NMER, NINE, LBDA, AAA ) 

CALL SPLINTt NMER ID, RAD I US, NSP, NMER, NINE, RAD, AAA) 

CALL SPL I NT (NMER ID, CUR V , NS P , NM ER , N I NE , CR V, AA A ) 

CALL SPLINTt NMERID, ALPHA , NS P , NMER , N I NE , AL PH , AAA ) 

CALL SPLINT (NMERID, BETA , NS P , NMER , N I NE , BT A , AA A ) 

110 DO 120 HM T= 1 , 3 

CPTIP(HMT) = 2,*CP*TIP(HMT) 

TWLMR = 2.*0MEGA*LBDA(4*HMT-3)-(0MEGA*RAD(4*HMT-3))**2 
TTIP = 1*-TWLMR/CPTIP(HMT) 
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120 WCR(HMT) = SQRT(TGROG*TTIP*TIP(HMT) ) 

DELMAX = WCR ( 2 ) / 20 • 

C CALCULATE COEFFICIENTS FOR MERIDIONAL VELOCITY GRADIENT EQUATION 
DO 125 1=1,9 

ALPH(I) = ALPH( I )/57. 295779 
CAL = COS ( ALPH (II) 

BTA(I) = BTA( I )/57. 295779 
SBETA = S I N l BT A ( I ) ) 

CBETA = COS ( BT A ( I ) ) 

A A ( I ) = CRVC I ) *CBETA**2— SBET A**2/RAD ( I )*CAL 
BBU) = -2«* OMEGA* SBET A* CAL 

ccm = CP* ( T IP ( I+1)-TIP(I ) ) -OMEGA* (LBDAf 1 4*1 )-LBDA III) 

125 CONTINUE 

C CALCULATE COEFFICIENTS FOR BLADE TO BLADE VELOCITY GRADIENT EQUATION 
DO 140 HMT= 1 • 3 
SAL = SIN(ALPH( 4*HMT - 3 ) ) 

TEMP = S A L*S I N ( BT A ( 4*HMT-3 ) ) /RAD ( 4* HMT- 3 ) 

EBB(HMT) = 2 • *OMEGA* S AL 
IFUX.E0.2) GO TO 132 
DO 130 BT B= 1 , 9 

130 A AB ( BT 8, HMT ) = TEMP + C S < HMT ) + ( CPR ( HMT ) -C S ( HMT ) ) *F LOAT ( BTB-1 ) /8 • 

GO TO 140 

132 DO 134 BT B= 1 , 9 

RC = l«/CS(HMT) + ( 1 • / CPR ( HMT)— 1»/CS( HMT ) ) * FLOAT ( BTB- 1 ) / 8* 

134 A AB ( BT B , HMT ) = TEMP+l./RC 
140 CONTINUE 

WEST = WTFL/RHOI P/ NM E R I D { NSP ) / ( NBB ( 1 ) +4 . *NBB < 2 ) +NBB ( 3 ) )*6. 

IF r WE ST.GT.WCR ( 2) ) WEST = WCR(2) 

WTFLSV = WTFL 

I F ( JZ* EQ. 3 ) WTFL = W C R ( 2 ) *RHO I P *NMER I D < NSP ) * ( NBB ( 1 ) +NBB ( 2 ) +NBB ( 3 ) ) 
IFIJZ.EQ.3) WEST = WCR<2) 

NCOUNT = 0 
145 IND = 1 
150 W ( 5 ) = WEST 
IP = 5 
IM = 5 

C CALCULATE W ON MID STREAMSHEET FROM MEAN TO HUB AND FROM MEAN TO TIP 
DO 160 1=1,4 
IP = IP+1 

WAS = W( IP-1 ) + <W( IP-1 )*AA( I P-1 )+BB( I P-1 ) )*DELTA+CC ( I P- 1 ) / W C I P-1 ) 
WASS = W( IP-1)+(WAS*AA{ I P ) +BB ( IP) )*DELTA+CC( IP-1)/ WAS 
W(IP) = ( W AS + W ASS ) /2 • 

IM = I M— 1 

WAS = W( IM+1 )-!W( IM+1 )*AAIIM+1 )+BB( IM+1) )*DELTA-CC ( IM)/W( IM+1) 

WASS = Wl IM+1)-(WAS*AA( I M) +BB ( IM) )*DELT A-CC ( I M ) / WAS 
160 W ( I M ) = < WAS+WASS >/2. 

CALCULATE RHO*W THROUGHOUT PASSAGE CROSS SECTION 

DO 180 HMT=1,3 
WBT B ( 5 ) = W(4*HMT— 3) 

WFIN(2,HMT) = WBTB ( 5 ) 

WSQ = W BT B ( 5 ) **2 

TWLMR = 2 • *QMEGA*LBDA ( 4*HMT-3 )— ( OMEGA*RAD( 4*HMT-3 ) ) **2 
TTIP = 1 • — ( WSQ+TWLMR )/CPTIPC HMT ) 

IF (TTIP.LT.O. ) GO TO 190 

RHOW = RHO I P ( HMT ) * ( 1 •— PLOSS ( HMT ) ) *TT I P**EXPON* WBTB ( 5 ) 

SUM ( HMT ) = RHOW 
DELTAB = NBB ( HMT ) / 8* 
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IP = 5 
IM = 5 

C CALCULATE W FROM BLADE TO BLADE ( HUB, MEAN, OR TIP) 

DO 17D 1=1,4 
IP = IP + 1 

WAS = WBTBI IP-1 > + ( WBTBI IP-1 ) *AAB( IP-1, HMT) +BBB I HMT ) )»DELTAB 
WASS = WBTBI IP-1 )+(WAS*AAB( IP, HMT )+BBBIHMT ) )*DELTAB 
WBTB(IP) = ( WAS+WASS ) /2. 

IM = IM-1 

WAS = WBTBI IM+1 )-( WBTBI IM+1 ) *AAB( I M+ 1 ,HMT ) +BBB l HMT )) «DEL TAB 
WASS = WBTBI IM+1 )-(WAS*AAB I I M , HMT ) +BBB I HMT ) )*DELTAB 
WBTBIIM) = I WAS+WASS ) /2. 

IE = IP 

165 WSQ = WBTBI I E ) **2 

TTIP = 1.-IWSQ+TWLMR )/CPTIP(HMT) 

IFITTIP.LT. 0. ) GO TO 190 

RHOW = RH0IP(HMT)*(1.-PL0SSIHMT) )*TTIP**EXPON*WBTBI IE) 

SUM (HMT) = SUM ( HMT )+RHOW 

IFII.EQ.4) SUM ( HMT ) = SUM I HM T ) -RHOW/ 2. 

IF! IE.EQ. IM) GO TO 170 
IE = IM 
GO TO 165 
170 CONTINUE 

WF I N 1 1 , HMT ) = WBTB(l) 

WF I N I 3 » HMT ) = WBTBI9) 

180 CONTINUE 

WTFLES = !SUM(1)*NBB( 1 ) +4. *SUM ( 2 ) *NBB I 2 ) +SUM I 3 ) *NBB ( 3 ) )/48.* 
1 NMERI9) 

IF! IND.GE.6. AND. ABSI WTFL-WTFLES) . LE . WTFL/ 1 0. E5 ) GO TO 200 
CALL CONT INI WE ST, WTFLES, IND , JZ , WTFL , DELMAX ) 

I F I IND.LT. 10) GO TO 150 
IF! IND.NE.10) GO TO 195 
WRI TE(6,2000) WTFLES 
GO TO 200 

190 WEST = WEST-DELMAX+3. 

NCOUNT = NCOUNT+1 
IFINCOUNT.LT. 50) GO TO 145 
195 WR I TEI 5 , 20 10 ) 

GO TO 240 

C SOLUTION HAS BEEN OBTAINED - PRINT SOLUTION 
200 DO 210 HMT =1,3 
DO 210 BT B= 1 » 3 

210 Vi WCR I B TB , HMT ) = WF IN I BTB ,HMT ) / WCR ( HMT ) 

WR I TE I 6 , 2020 ) 

WRI TEI 6,2040) I I WF INI BTB, HMT > , BTB= 1 , 3 ) , HMT = 1 , 3 ) 

WRITE! 6,2030) 

WRITE! 5,2040) I tWWCRI BTB,HMT),BTB=1,3),HMT=1,3) 

C CALCULATE MASS FLOW DISTRIBUTION AND STREAMLINE SPACING 
DO 220 1=1,3 

Y I I ) = SUM! I ) *NB6 I IJ/8.0 
220 XII) = FLOAT ( 1-1 )/2. 

CALL PABC (X, Y.A.B.C) 

SL190 = N MER I 9 ) / 100. 

YAVE = WTFLES/NMERI9 ) 

WRITEI6.2050) 

DO 230 1=1,11 

FRSLD = FLOAT! I-D/10. 

WFFK = A*FRSLD**3/3. +B*FRSLD**2/2.+C*FRSLD 
WFFR = WFFR/YAVE 
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RSLSP = YAVE/ ( A*FRSLD**2+B*FRSLD+C ) 

8100 = RSLSP*SL100 

230 WR I TEI 5 » 206U ) FRSLD, WFFR , RSLSP, 8100 
240 WRITE(6,1000) 

WTFL = WTFLSV 
GO TO 10 

1000 FORMAT ( 1H1 ) 

1010 FORMAT ( 1X»80H 
1 

1020 FORMAT ( 1615) 

1030 FORMAT l 8F 10 . 5 ) 

1040 FORMAT < 1X.8G16.5) 

1100 FORMAT ( 25H JX JZ KR1 KR2 NSP) 

1110 FORMAT (8 X.3HGAM, 14X, 2HAR , 1 IX , 5H0MEGA , 12X , 4HWTFL ) 

1120 FORMAT ( 1HL , 1 IX , 5HRH0 I P , 1 IX , 5HPL0SS , 12X , 3HNBB , 14X , 2HCS , 1 3X , 3HC PR ) 
1130 FORMAT ( 6H HUB .5G16.5/6H MEAN .5G16.5/6H TIP ,5G16.5) 

1140 FORMAT ( 1HL,9X» 12HNMERID ARRAY) 

1150 FORMAT { 10X, 12HTEMPIP ARRAY) 

1160 FORMAT ( 10X»12HLAMBDA ARRAY) 

1170 FORMAT ( 10X, 1 2HRAD I US ARRAY) 

1180 FORMAT ( 10X, 10HCURV ARRAY) 

1190 FORMAT ( 10X, 1 1HALPHA ARRAY) 

1200 FORMAT t 10X, 10HBETA ARRAY) 

1500 FORMAT ( 39HLNSP MUST BE GREATER THAN OR EQUAL TO 2) 

1510 FORMAT (30HLNSP CANNOT BE GREATER THAN 25) 

1520 FORMAT ( 26HL JX MUST BE EITHER 1 OR 2) 

1530 FORMAT ( 49HL WHEN JX = 2, CS ANO CPR MUST HAVE THE SAME SIGN) 

2000 FORMAT ( 42HLTHE PASSAGE IS CHOKED WITH A MASS FLOW OF.G16.5) 

2010 FORMAT ( 46HLA SOLUTION CANNOT BE OBTAINED AT THIS STATION) 

2020 FORMAT ( 1HL* 10X, 2 5HVE LOCITIES ACROSS PASSAGE) 

2030 FORMAT! 1HL , 10X , 39HCR I T I C AL VELOCITY RATIOS ACROSS PASSAGE) 

2040 FORMAT (9X,7HSUCT ION, 1 IX, 3HMID, 13X, 8HPRESSURE/6H HUB ,3G16.5/ 

1 6H MEAN ,3G16.5/6H TIP .3G16.5) 

2060 FORMAT ( 3X» 5G25.5) 

2050 FORMAT { 105HL FRACTIONAL STREAMLINE MASS FLOW 

1 RATIO OF ACTUAL TO AVERAGE NORMAL THICKNESS B /15X, 

2 BHD I STANCE, 17X,8HFRACTION, 13X, 18HSTREAMLINE SPACING, 6X, 

3 2 1H { FOR 100 STREAMLINES) ) 

END 


Dictionary of Variables in Main Program 

o 

A coefficient of X calculated by subroutine PABC 

AA array of coefficients a n from eq. (A2) 

AAA array of dummy values not used in calculation 

AAB array of coefficients a^ from eq. (A4) 

ALPH array of values of a at eight equal intervals from hub to tip 

ALPHA input array 

AR input 
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B 

B100 

BB 

BBB 

BETA 

BTA 

BTB 

C 

CAL 

CBETA 

CC 

CP 

CPR 

CPTIP 

CRV 

CS 

CURV 

D ELM AX 

DELTA 

DELTAB 

EX PON 

FRSLD 

GAM 

HMT 

I 

IE 

IM 

IND 

IP 

JX 


coefficient of X calculated by subroutine PABC 
normal streamline spacing for 100 streamlines 
array of coefficients b n from eq. (A2) 
array of coefficients b^ from eq. (A4) 
input array 

array of values of /3 at eight equal intervals from hub to tip 
integer subscript used to indicate blade -to -blade position 
constant coefficient calculated by subroutine PABC 
cos ot 
cos /3 

array of coefficients c fi from eq. (A2) 
specific heat at constant pressure, c 

Sr 

input array 

array of values of 2CpTJ 

array of values of l/r c at eight equal intervals from hub to tip 
input array 
input array 

maximum change in WEST permitted in subroutine CONTIN 

one -eighth the distance from hub to tip 

one -eighth the blade -to -blade streamline normal distance 

i/(y-D 

fractional distance from hub to tip 
input 

integer subscript used to indicate hub, mean, or tip position 
DO loop index 
integer subscript 
integer subscript 

integer switch - used for logical control in subroutine CONTIN 

integer subscript 

input 
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JZ 

KR1 

KR2 

LAMBDA 

LBDA 

NBB 

N COUNT 

NINE 

NMER 

NMERID 

NSP 

OMEGA 

PLOSS 

RAD 

RADIUS 

RC 

RHOIP 

RHOW 

RSLSP 

SAL 

SBETA 

SL100 

SRW 


SUM 

TEMP 

TEMPIP 

TGROG 


input 

input 

input 

input array 

array of values of X at eight equal intervals from hub to tip 
input array 

number of times that a solution could not be obtained because of excessive 
velocities 

9 

array of distances from hub 

input array 

input 

input 

input array 

array of values of r at eight equal intervals from hub to tip 
input array 

( r c)b-b 
input array 

PW 

ratio of actual to average streamline spacing 
sin a 
sin /3 

one -hundredth the distance from hub to tip 

integer code variable that will cause SPLINT to write out data useful for de- 
bugging (If SRW = 16, SPLINT will write input and output data. SRW is not 
an input item . ) 

array of integrated values of PW at hub, mean, and tip 
temporary storage 
input array 
2yR/(y+l) 
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TIP 


array of values of T? at eight equal intervals from hub to tip 
TTIP T/TJ 

TWLMR 2wX - (cor) 2 

W array of values of W ^ at eight equal intervals from hub to tip 

WAS initial estimate for W at next point 

WASS secondary estimate for W at next point 

WBTB array of values of W at eight equal intervals from blade to blade 

WCR array of values of W cr at hub, mean, and tip 

WEST estimated value for W ^ „ 

mid, mean 

WFFR weight flow fraction between some point and hub 

WFIN array of final calculated velocities to be printed as output 

WSQ W 2 

WTFL input 

WTFLES calculated mass flow if W . , „ = WEST 

mid, mean 

WTFLSV used to save input value of WTFL when choking mass flow is to be calculated 
WWCR array of values of W/W cr for output 

X array of fractional distances from hub 

Y array of values of mass flow per unit radial distance 

YAVE average value of mass flow per unit radial distance 


Subroutine CONTIN 


Subroutine CONTIN finds the value of W^^ mean which will give either the desired 


value of mass flow WTFL or the choking mass flow, depending on the input value of JZ. 

The actual calculation of mass flow corresponding to a given value of W m ^ mean is 

done by the main program ; CONTIN determines the value of W ,, which should 

iii iq ^ in G3.n 

be used. If WTFL is larger than the choking mass flow, the choking mass flow will be 
automatically calculated, and a message to this effect will be printed. 

The input arguments for CONTIN are as follows: 

XEST estimated value of W . , 


x " mid, mean 

YCALC calculated mass flow for W . j 


mid, mean 


= XEST 
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IND 


switch controlling action to be taken by CONTIN 
JZ program input, see fig. 2 

YGIV desired mass flow, WTFL, when JZ = 1 or 2 
XDEL maximum change in XEST in one iteration 
The output arguments for CONTIN as follows: 

XEST estimated value of W m ^ mean ^ or nex * iteration 

IND switch indicating whether a solution can or cannot be obtained by CONTIN 

The internal variables in CONTIN are as follows: 

2 

APA coefficient of X calculated by subroutine PABC 

BPB coefficient of X calculated by subroutine PABC 

CPC constant coefficient calculated by subroutine PABC 

DISCR discriminant, |/b^ - 4ac, of parabola obtained by subroutine PABC 

NCALL number of times CONTIN has been called for a given case 

X array of three values of XEST from previous calls 

Y array of three values of YCALC from previous calls 


SUBROJT INE CONT INI XEST, YCALC, I NO , J Z , YG I V , XDE L ) 
DIMENSION X(3),Y(3) 

NCALL = NCALL+1 

I F ( IND.NE. I. AND. NCALL. GT. 50) GO TO 160 
GO TO ( 10,30,40,50,60,80,130), IND 
C FIRST CALL 
10 NCALL = 1 

IFIYCALC. GT. YGIV. AND. JZ.EQ.l) GO TO 20 
IND = 2 
Y ( 1 ) = YCALC 
X ( 1 ) = XEST 
XEST = XEST+XDEL 
RETURN 
20 IND = 3 

Y ( 3 ) = YCALC 
X ( 3 ) = XEST 
XEST = XEST-XDEL 
RETURN 

C SECOND CALL 
30 IND = 4 

Y l 2 ) = YCALC 
X ( 2 ) = XEST 
XEST = XEST+XDEL 
RETURN 
40 IND = 5 

Y ( 2 ) = YCALC 
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X ( 2 ) = XE ST 
XEST = XEST-XDEL 
RETURN 

C THIRD OR LATER CALL - FIND SUBSONIC OR SUPERSONIC SOLUTION 
50 Y ( 3 ) = YCALC 
X ( 3 ) = XEST 
GO TO 70 
60 Y ( 1 ) = YCALC 
X ( 1 ) = XEST 

70 IF( YGIV.LT.AMINKYt 1 ) ,Y(2)*Y(3>)) GO TO (90*95), JZ 
75 IND = 6 

CALL PABC (X, Y, APA,BPB,CPC) 

D I SCR = BPB**2-4.*APA*(CPC-YGIV) 

I F ( DISCR.LT.O. ) GO TO 110 

XEST = -8 P B-S I GN ( SQR T ( D I SCR ) * APA ) 

I F ( JZ.EQ.2.AND.APA.LT.0. ) XEST = -BPB-SQRT ( D I SCR ) 

XEST = XEST/ 2 • /APA 
IF(XEST.GT.X(3)+XDEL) GO TO 95 
IFIXEST.LT. X(l)-XDEL) GO TO 90 
RETURN 

C FOURTH OR LATER CALL - (NOT CHOKED) 

80 IF( XEST.GT.X(3) ) GO TO 50 
IF( XEST.LT.X( 1 ) ) GO TO 60 

Y ( 2 ) = YCALC 
X ( 2 ) = XEST 
GO TO 70 

THIRD OR LATER CALL - SOLUTION EXISTS? 

BUT RIGHT OR LEFT SHIFT REQUIRED 
90 IND = 5 
C LEFT SHIFT 

Y ( 3 ) = Y ( 2 ) 

X ( 3 ) = X ( 2 ) 

Y ( 2 ) = Y ( 1 ) 

X ( 2 ) = X( 1) 

XEST = X( D-XDEL 
RETURN 

95 IND = 

C RIGHT SHIFT 

Y ( 1 ) = Y ( 2 ) 

X( 1) = X ( 2 ) 

Y ( 2 ) = Y ( 3 ) 

X ( 2 ) = X ( 3 ) 

XEST = X ( 3 ) +XDEL 
RETURN 

C THIRD OR LATER CALL - APPEARS TO BE CHOKED 
110 XEST = -BPB/2./APA 
IND = 7 

IF (XU) . LE. XEST. AND. XEST.LE.Xt 3) ) RETURN 
IF(XEST.LT.Xd) ) GO TO 90 
GO TO 95 

C FOURTH OR LATER CALL - PROBABLY CHOKED 
130 IF ( YCALC. GE.YGIV ) GO TO 80 
IND = 10 
RETURN 

C NO SOLUTION FOUND IN 50 ITERATIONS 
160 IND = 11 
RETURN 
END 
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Subroutine SPLINT 


This subroutine calculates interpolated values using a cubic spline curve fitted 
through the given data. The input variables for SPLINT are as follows: 

X array of spline point ordinates 

Y array of function values at spline points 

N number of X and Y values given 

Z array of ordinates at which interpolated function values are desired 

MAX number of Z values given 
The output variable for SPLINT is as follows: 

YINT array of interpolated function values 

If SRW = 16 in COMMON, or if some element of the z -array falls outside of the 
interval for the x -array, input and output data for SPLINT are printed. This is useful 
in debugging. 


SUBROUTINE SPLINT ( X , Y , N , Z , M AX , Y I NT , 0 YDX ) 

C 

C SPLINT CALCULATES INTERPOLATED POINTS AND DERIVATIVES 
C FOR A SPLINE CURVE 

C END CONDITION - SECOND DERIVATIVE AT EITHER END POINT IS ONE-HALF 
C THAT AT THE ADJACENT POINT 
C 

COMMON SRW 

DIMENSION X(N) f Y (N ), Z ( MAX) , Y INTI MAX) , DYDX(MAX) 

DIMENSION G( 100) ,SB( 100) »EM( 100) 

INTEGER SRW 
IF(MAX.LE.O) RETURN 
III = SRW 
SB ( 1 ) = -.5 
G { 1 ) = 0 
N0=N-1 

IFINO.LT. 2) GO TO 20 
DO 10 1=2, NO 
A = ( X ( I)-X( I-l) )/6. 

C = <X( I+l)-X( I ) )/6. 

W = 2 . * ( A + C ) -A*SB ( I-l ) 

SB ( I ) = C/W 

F = ( Y ( 1 + 1 )-Y( I ) )/ (X< 1 + 1 )-X( I ) )-<Y( I )-Y( I-l ) )/(X( I )-X(I-l ) ) 

10 G ( I ) = ( F — A*G ( I-l) )/W 

20 EMIN) = GIN— l)/(2.+SB(N-l) ) 

DO 30 1=2, N 
K = N+l-I 

30 EM ( K ) = G(K)-SBtK)*EM(K+l) 

DO 140 1=1, MAX 
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K = 2 

I F { Z ( I ) —X ( 1 ) ) 70,60,90 
60 Y I N F ( I )=Y( 1) 

SK = X(K)-X(K-1) 

GO TO 130 

70 1F{Z(I).GEJ1.1»X(1)-.1*X(2))) GO TO 120 
WRITE (6,1000) Z(I) 

SR W = 16 
GO TO 120 
80 K=N 

IF{Z(II.LE,lltl*X{N)-.l*X(N-l)n GO TO 120 
WRITE (6, 1000 ) Z ( I ) 

SR W = 16 
GO TO 120 

90 IF(Z(I )-X(K) ) 120, 100,110 

100 Y I N T ( I ) = Y ( K ) 

SK - X(K)-X(K-1) 

GO TO 130 
110 K=K + 1 

IF(K-M) 90,90,80 
120 CONTINUE 

SK = X(K)-X(K-1) 

YINT(I) = EM(K-1)*(X(K)-Z( I) )**3/6./SK +EM ( K * ( Z ( I i - X ( K- 1 ) ) **3/ 6 . 

1 / SK+ ( Y ( K ) / SK -EM(K)*SK / 6 . ) * ( Z ( I ) -X ( K- 1 ) ) + ( Y ( K- 1 ) / SK -EM(K-l) 

2 * SK/ 6 • )*(X(K)-Z(I) ) 

130 D YOX ( I )=-EM(K-l)*(X(K)-Z (I i ) **2/ 2.0/SK +EM ( K ) * (X(K-l)-Z ( I > >**2/2. 

1 /SKM Y< K)-Y (K-l ) ) / SK - ( EM ( K ) -EM { K-l ) )*$K/6. 

140 CONTINUE 

MX A = M AXO ( N , MAX ) 

IF ( SRW. EO. 16) WRITE( 6, 1010 ) N, MAX • (X(I),Y(I),Z(I),Y[NT(I) ,DYDX(I ) , 
1 1=1, MXA) 

SRW =111 
RETURN 

1000 FORMAT (54H SPLINT USED FOR EXTRAPOLATION. EXTRAPOLATED VALUE = , 
1G14.6) 

1010 FORMAT ( 2 X , 2 1 HNO • OF POINTS GIVEN =,I3,30H, NO. OF INTERPOLATED PO 
1 I NTS =, I3/10X, 1HX, 19X, 1HY, 16X, 1 1HX- INTERPOL 9X, 1 1HY- I NTERPOL • , 

28 X, 14HD YD X- INTERPOL. / ( 5G20.8 ) ) 

END 


Subroutine PABC 

p 

Subroutine PABC calculates coefficients A, B, and C of the parabola y = Ax + 
Bx + C passing through three given X, Y points. 
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SUBROUTINE P ABC ( X , Y, A , B, C > 

C SUBROUTINE P ABC CALCULATES COEFFICIENTS A, B, C OF THE PARABOLA 
C Y = A*X**2 + B*X + C, PASSING THROUGH THREE GIVEN X,Y POINTS 
DIMENSION X { 3 > ,Y(3) 

Cl = X( 3 » -X C 1 ) 

C2 = (Y(2)-Y( I) )/(X( 2 ) — X ( 1 ) ) 

A = ( C 1 *C2— Y (3I+YIL) )/CI/(X(2)-X(3)I 

B = C2 - ( X I 1 ) +X t 2 ) >*A 

C = Y( 1 ) — X ( 1 )*B— X{ 1 ) * * 2 * A 

RETURN 

END 


Lewis Research Center, 

National Aeronautics and Space Administration, 
Cleveland, Ohio, October 27, 1970, 

720-03. 



APPENDIX A 


ANALYTICAL EQUATIONS AND SOLUTION PROCEDURE 

Two ordinary differential equations are utilized to determine the velocity gradients 
in both the blade -to -blade and the hub-to-tip directions. These equations are both de- 
rived in appendix B from a general velocity gradient equation given in reference 5. The 
hub-to-tip velocity gradient equation is 


dW 

dn 


aw + b 

n n 



(Al) 


where 


'\ 


dn 


dn 


2 2 
cos (3 _ cos a sin /3 

a n 


b fi = -2 co cos cl sin /3 f 


dT? .. 
c „ " P — -■»- 


J 


(A2) 


An analytical solution of equation (Al) cannot be obtained, since /3 is not given as 

an analytical function. However, a numerical solution can be readily obtained. The 

numerical method used is the Heun method (ref. 7, p. 236). 

With an estimated value for mean equation (Al) can be integrated from 

mean to hub and from mean to tip by the Heun method. This gives the entire midchannel 

velocity W m ^ distribution from hub to tip at the given station for the estimated value 

of W .. „ . The first estimate for W ., depends on the input value of JZ. 

mid, mean mid, mean r * 

If JZ = 1 or 2, W . . is estimated based on one -dimensional flow using the channel 

miQj in 

cross section area, input mass flow, and inlet stagnation density p!. If JZ = 3, 

W .j is estimated to be W „ at the mean between hub and tip. 

mid, mean cr 

The blade -to -blade velocity gradient equation (derived in appendix B) is 


dW 


dn. 


= a b W + b b 


(A3) 


b-b 
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where 


*b = 


sin a sin ff 
r 


1 

( r c>b-b ► 


b^ = 2lo sin a 


(A4) 


Equation (A3) could be solved analytically if (r or l/(r c ) bb is assumed to 

vary linearly along the blade -to -blade orthogonal. However, this solution is very com- 
plex and involves an integral which must be evaluated numerically. Therefore, a numer- 
ical solution of equation (A3) requires much less computation and will give at least as 
good accruacy as the analytical solution. The numerical method used in the program is 
the Heun method (ref. 7, p. 236). The solution depends on whether JX = l(l/(r c ) b b is 
assumed to vary linearly from blade to blade) or JX = 2 (( r c )b_k is assumed to vary 
linearly from blade to blade). The value of W ^ to be used is obtained from the nu- 
merical solution to equation (Al). Then equation (A3) is integrated numerically from the 
midchannel to each blade surface separately at hub, mean, and tip. 

With the velocities obtained, the weight flow past the channel cross section can be 
computed from 



PW dn b _ b 


dn 


(A 5) 


integrating from blade to blade and from hub to tip. 

The density p is calculated using the equation for isentropic flow, with a correction 
for loss in the relative stagnation pressure (given as PLOSS in input). This equation for 
density is 


P = P^l - 


W 2 + 2 to A - (qjr) 2 

2c p T i 


i l/y-i 


l - 


p” 
los s 

p." 
lseny 


(A6) 


(This equation is derived in appendix B.) 

The inner integrals are computed at hub, mean, and tip using trapezoidal integration 
over eight equal intervals from suction to pressure surface. 

After the inner integral is computed in equation (A5) at hub, mean, and tip, the outer 
integral is approximated by Simpson's rule. This results in a value for the mass flow w 


based on the estimated value for W 


mid, mean' 


This value of w can be compared with 
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the input value WTFL. The procedure for calculating the mass flow for a given value of 
W mi d m ean is outlined graphically in figure 8. 

If JZ = 1 or 2, further calculations are made to determine a value of W m ^ mean 
that will give a mass flow w close to the input value WTFL. If WTFL is less than the 
choking mass flow, there are two values of W ^ mean which will give w = WTFL. If 
JZ = 1, the smaller, or subsonic, solution will be found. If JZ =2, the larger, or super- 
sonic, solution will be found. When |w - WTFL| < 10" *WTFL the calculations are 
stopped, and the solution is printed. If WTFL is larger than the choking mass flow, no 
solution exists. In this case, the choking mass flow is calculated, and this solution is 
printed with the message "THE PASSAGE IS CHOKED WITH A MASS FLOW OF XXXX". 

If JZ = 3, further calculations are made to determine the choking mass flow, regard- 
less of the input value of WTFL. The second and third estimates for W mid mean are 
each increased from the previous estimate by 5 percent of W cr at the mean (the first 


estimate was W„ 


are calculated to 


• =W1. New estimates for W . , 

mid, mean cr' mid, mean 

obtain the maximum weight flow. These new estimates are based on three previous ap- 
proximations for W m ^ d mean with the corresponding mass flows. When the new esti- 
mate for W m j, d mean is estimated to be within the range of the three values used for 
the last quadratic fit, the corresponding mass flow is calculated and this value is ac- 
cepted as the choking mass flow. 

The final output for a case gives information concerning the mass flow distribution 
along the meridional normal. A quadratic fit is used from hub to tip (corresponding to 
Simpson's rule in integrating eq. (A5) from hub to tip). The quadratic fit determines 
the mass flow fraction at 10 equally spaced intervals from hub to tip, as well as the 
streamline spacing. 
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APPENDIX B 


DERIVATION OF EQUATIONS 


A general velocity gradient equation can be derived from the force equation. If it is 
assumed that entropy is constant along an arbitrary space curve, with the distance along 
the curve denoted by q, the following equation can be obtained: 


dW 

dq 


a^ + b^ + 

dq dq 




Mb' 


dq W \ dq 



(Bl) 


where 


as Wcos tt cos 2 g _^nJ^ slnacosg ^_ 2 ^ 3tng 
r„ r dm 


, W sin a cos 2 P R dW m 

b = - + cos a cos p 

r dm 


r ( B2 ) 


c = W sin a sin P cos P + r cos 


0 


, dm 


+ 2 u) sin c* 


Equations (Bl) and (B2) are the same as equations (B13) and (B14) of reference 5. 

To derive equations (Al) and (A2), it is assumed that q is equal to n, the distance 
along an orthogonal to the streamlines in a meridional plane. Then q = n, dr/dq = dr/dn 
cos oi, dz/dq = dz/dn = -sin ot, and d0/dq = 0. Since it is assumed that c^ is constant, 
equations (Bl) and (B2) reduce to equations (Al) and (A2). 

To derive equations (A3) and (A4), it is assumed that q is equal to n b b , the dis- 
tance along a blade -to -blade streamline orthogonal. Along the blade to blade normal, 
then, r and z are both functions of n b _ b - Alternately, r and z can be considered to 
be functions of m, where m is a function of n^ ^. Then we can obtain the following 
equations : 


dr _ dr dm 
dn b _ b dmdn b _ b 


- sin a sin P 
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dz _ dz dm 
dn b-b dm dn b-b 


- cos oi sin P 


Also, 


r — — — = cos P 


dn 


b-b 


In the blade -to -blade direction uniform upstream conditions are assumed; therefore, 


dh! 


dA 


= 0 


dn b-b dn b-b 


With these substitutions and simplification, equations (Bl) and (B2) become 


dW dW 

d ^ _ W sin a sin P + s j n a _ s ^ n p cos p + cog 2 p £ (B3) 


dn 


b-b 


dm 


dm 


Further simplification of equation (B3) is possible by evaluating the derivatives 
dW^/dm and dW 0 /dm in terms of angles and blade -to -blade streamline curvatures. 


m' 


We make use of the fact that W 0 = W sin P and W m = W cos P (see fig. 3). Also 
dj3/dm = (d|3/ds)/cos P and d|3/ds is the blade -to -blade streamline curvature 1/ ( r c )| 3 _) ;) * 
This gives 


'A 


dWg 

dm 


dW 


m 


dm 


= cos P 


dW 

W 

dm 

( r c^b-b 

dW _ 

W tan P 

dm 

< r o>b-b 






When equations (B4) are substituted in equation (B3) we obtain 


dW _ W sin a sin P 


W 


dn- 


b-b 


(rJ 


. + 2 a) sin a 


c'b-b 


(B4) 


(B5) 
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Equations (A3) and (A4) are obtained directly from this equation. 

The density corresponding to a given velocity may be calculated using the isentropic 
relation with a correction for loss in total pressure. The equation for this is equa- 
tion (A6), which is derived here. We start with 




and use each of the following relations: 


- = (-) 

P" VT”/ 




p» = _E1 


RT' 


T N l/(r-l) / T \l/(y-i) /Tiy/(r-l) 

rp I f 




T i> 


t rn T v 1/ (y 1) Qt 

— = — 
T"/ 


p." 

isen 


Using these relations results in 


p.” RT” = p." 
isen ^isen 


n" = n." - p ,T 

1 *isen *Toss 




i/(y-i) 


l - 


loss 


isen; 


(B6) 


From equation (3) of reference 5 


T_ = j _ W 2 + 2coA - (cur) 2 


T i 


2o P T i 


When this is substituted in equation (B6), equation (A6) is obtained. 
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JX 

JZ 

/ 

/ 

GAM 

A 4 


1 5|6 10 1 II 15 |l6 20 121 25| 26 3o|31 


TITLE (first case only) 




LaBEL (every case) 


RHOIP 


.00230SS7 


RHOIP 


.00030567 


KR1 

A 

AR 

[Z7/6.67 
Hub 

PLOSS 

| .04*31 
Mean 

PLOSS 
.0482 


KR2 


NSP 


JX 


3 

OMEGA 

439- 

NBB 

. 042 09 

NBB 

.0545 


4o| 41 


fl -Lii 
! [2-Lii 


Linear curvature 
Linear radius of curvature 

WTFL 
. 0/839 


cs 

- 55 : 7 / 

cs 

l /4. /<o 


50 1 51 


} KR1= {l- 

CPR 

-8.QZ 

| CPR 
l -8.42 


Omit this card 
Supply this card 


JZ = 


1 - Subsonic solution 

2 - Supersonic solution 

3 - Choked flow solution 


KR2 = 


0 - Omit these cards 

1 - Supply one card 

3 - Supply three cards 


Tip 


RHOIP 

PLOSS 

NBB 

CS 

CPR 


.00230557 

■ 0481 

.0715 

-22.0 

-6.30 



Use NAMELIST input for NMERID, TEMPIP, LAMBDA, RADIUS, CURV, ALPHA and BETA 


JNAMi _/V/agRID- 0.^875,. 37^ TfMPIP = 3* 518. U, /AMBDA = 3 * 1024 ., PADS US '.875, /. 06 zr, 
BETA x - 76.3, -zes _ A 


1.25, 


Figure 2. - CHANEL input form. 



Figure 3. - Cylindrical coordinate system and velocity components. 
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Reference 
0 coordinate 


m 



Figure 6. - Blade-to-blade channel layouts. Angular stacking points, (T) , © , and @ are located from figure 5. Input items BETA, NBB, CS, and CPR must be given as true angle, length, and 
curvatures, not as quantities measured in m,0 plane. Correct quantities can be calculated from measurements in m, 0 plane. 



SAMPLE PROBLEM 
NEGATIVE REACTION TURBINE 


JX 

JZ KR1 

KR2 

NSP 



1 

1 1 

3 

3 




GAM 


AR 

OMEGA 

WTFL 


1.40000 


1716.57 

488.000 

0.1 8 390E-01 



RHOI P 

PLOSS 

NBB 

cs 

C PR 

HUB 

0. 23056t-02 

0 .482006—0 1 

0.42080E-01 

-15.7100 

-8.82000 

MEAN 

0. 23056E-02 

0.48200E-01 

0. 54500E-01 

-14.1600 

-8.42000 

TIP 

0.230566-02 

0 • 48 200E— 0 I 

0.71 500E-0 1 

-12.0000 

-6.30000 


NIMERIO ARRAY 
0 0.18750 

TEMP IP ARRAY 
518.670 518.670 

lambda array 

1024.00 1024.00 

RADIUS ARRAY 
0.87600 1.0625C 

C UR V ARRAY 
0 0 

ALPHA ARRAY 
0 0 

BETA ARRAY 

-5. 00300 -16.3000 


0. 37500 
518.670 
1024.00 
1.26000 
0 
U 

-28.5000 


HUB 
MEAN 
T I P 


VELOCITIES 

SUCTION 

887.531 

978.047 

1085.45 


ACROSS PASSAGE 
M I U 

662.924 
691. 750 
744. 150 


PRESSURE 
531. 104 
528.942 
564.698 


HUB 
MEAN 
T IP 


CRITICAL VELOCITY RATIOS ACROSS PASSAGE 


suction 
0. 73633 
1.02137 
1. 12304 


M ID 

0.69780 
0.72239 
0 . 76992 


PRESSURE 
0.55904 
0.55237 
0. 58426 


hRACTIONAL STREAMLINE 
DISTANCE 
0 

0. 10000 
O.2UO0O 
0. 30000 
0 . 4U0U0 
0. 5UOOO 
0.60000 
0. 70000 
0. 80U00 
0.90000 
1.00000 


MASS FLOW 

fraction 

0 

0. 72987E-01 
0.15043 
0.23292 
0.32102 
0.41532 
0.51639 
0.62481 
0.741 17 
0.86604 
1.00000 


RATIO OF ACIUAL TO AVERAGE NORMAL 
STREAMLINE SPACING (FOR 100 

1.40942 0 

1.33120 0 

1.25205 0 

i. 17376 0 

1.09766 0 

1.02472 0 

0.95555 0 

0.69051 0 

0.62976 0 

0.7/329 0 

0.72099 0 


Figure 7. - Output listing. 


THICKNESS B 
STRtAML INES ) 
52853E-C2 
49920E-02 
46952E-02 
4401 6E-02 
4 U62E-G2 
3 842 7E-G2 
3 58 3 3E-02 
3J396E-C2 
3111 6E-U2 
28998E-02 
2 703 7E-02 
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(1) Estimate mearr 

(2) Calculate W mjd from mean to hub and mean to tip from eq. (Al). 

(3) Calculate W from midchannel to suction and pressure blade surfaces at 
hub, mean, and tip using eq. (A3). 

(4) Calculate total mass flow for channel cross section using eq. (A5). 

(5) Repeat steps (1) to (4) with new estimates for W mid mean until either 
desired mass flow or choking mass flow is obtained.’ 

Figure 8. - Calculation procedure chart. 
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